#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define MAX_LINE_LEN 1024
#define EPS 1e-10
#define MAX_ITER 1000

// Read CSV file and extract numeric columns
double** read_csv_numeric(const char *filename, int *rows, int *cols) {
    FILE *fp = fopen(filename, "r");
    if (!fp) {
        printf("Failed to open file: %s\n", filename);
        printf("Please check if the file exists and the path is correct\n");
        return NULL;
    }

    char line[MAX_LINE_LEN];
    int col_count = 0;

    // Read header line to count columns
    if (fgets(line, MAX_LINE_LEN, fp)) {
        char *token = strtok(line, ",");
        while (token) {
            col_count++;
            token = strtok(NULL, ",");
        }
    }

    // Calculate numeric columns (skip first 3 columns)
    *cols = col_count - 3;
    if (*cols <= 0) {
        printf("No valid numeric columns found in the file\n");
        fclose(fp);
        return NULL;
    }

    // Count total rows
    *rows = 0;
    while (fgets(line, MAX_LINE_LEN, fp)) {
        (*rows)++;
    }

    // Allocate memory for data
    double **data = (double**)malloc(*rows * sizeof(double*));
    for (int i = 0; i < *rows; i++) {
        data[i] = (double*)malloc(*cols * sizeof(double));
    }

    // Read data content
    rewind(fp);
    fgets(line, MAX_LINE_LEN, fp); // Skip header line

    int row = 0;
    while (fgets(line, MAX_LINE_LEN, fp) && row < *rows) {
        int c = 0;
        char *token = strtok(line, ",");
        int numeric_col = 0;
        while (token) {
            // Skip first 3 columns (Area, Item, Year)
            if (c >= 3) {
                data[row][numeric_col++] = atof(token);
            }
            c++;
            token = strtok(NULL, ",");
        }
        row++;
    }

    fclose(fp);
    return data;
}

// Center data by subtracting mean
void center_data(double **data, int rows, int cols, double *mean) {
    // Calculate mean for each column
    for (int j = 0; j < cols; j++) {
        mean[j] = 0.0;
        for (int i = 0; i < rows; i++)
            mean[j] += data[i][j];
        mean[j] /= rows;
    }

    // Subtract mean from each data point
    for (int i = 0; i < rows; i++)
        for (int j = 0; j < cols; j++)
            data[i][j] -= mean[j];
}

// Calculate covariance matrix
void covariance_matrix(double **data, int rows, int cols, double **cov) {
    for (int i = 0; i < cols; i++) {
        for (int j = 0; j < cols; j++) {
            double sum = 0.0;
            for (int k = 0; k < rows; k++)
                sum += data[k][i] * data[k][j];
            cov[i][j] = sum / (rows - 1);
        }
    }
}

// Jacobi algorithm for eigenvalue decomposition
void jacobi_eigen(double **A, double *eigenvalues, double **eigenvectors, int n) {
    // Initialize eigenvalues and eigenvectors
    for (int i = 0; i < n; i++) {
        eigenvalues[i] = A[i][i];
        for (int j = 0; j < n; j++)
            eigenvectors[i][j] = (i == j) ? 1.0 : 0.0;
    }

    // Jacobi iterations
    for (int iter = 0; iter < MAX_ITER; iter++) {
        // Find the maximum off-diagonal element
        int p = 0, q = 1;
        double max = fabs(A[p][q]);
        for (int i = 0; i < n; i++) {
            for (int j = i + 1; j < n; j++) {
                if (fabs(A[i][j]) > max) {
                    max = fabs(A[i][j]);
                    p = i; q = j;
                }
            }
        }

        // Convergence check
        if (max < EPS) break;

        // Calculate rotation angle
        double theta = 0.5 * atan2(2 * A[p][q], A[q][q] - A[p][p]);
        double c = cos(theta), s = sin(theta);

        // Update matrix A
        double app = c*c*A[p][p] - 2*s*c*A[p][q] + s*s*A[q][q];
        double aqq = s*s*A[p][p] + 2*s*c*A[p][q] + c*c*A[q][q];
        A[p][p] = app;
        A[q][q] = aqq;
        A[p][q] = A[q][p] = 0.0;

        // Update other elements
        for (int j = 0; j < n; j++) {
            if (j != p && j != q) {
                double apj = c*A[p][j] - s*A[q][j];
                double aqj = s*A[p][j] + c*A[q][j];
                A[p][j] = A[j][p] = apj;
                A[q][j] = A[j][q] = aqj;
            }
        }

        // Update eigenvectors
        for (int j = 0; j < n; j++) {
            double vpj = c*eigenvectors[p][j] - s*eigenvectors[q][j];
            double vqj = s*eigenvectors[p][j] + c*eigenvectors[q][j];
            eigenvectors[p][j] = vpj;
            eigenvectors[q][j] = vqj;
        }
    }

    // Extract eigenvalues
    for (int i = 0; i < n; i++)
        eigenvalues[i] = A[i][i];
}

// Perform PCA and output results
void pca(double **data, int rows, int cols) {
    // Allocate memory
    double *mean = (double*)malloc(cols * sizeof(double));
    double **cov = (double**)malloc(cols * sizeof(double*));
    for (int i = 0; i < cols; i++) {
        cov[i] = (double*)malloc(cols * sizeof(double));
    }

    double *eigenvalues = (double*)malloc(cols * sizeof(double));
    double **eigenvectors = (double**)malloc(cols * sizeof(double*));
    for (int i = 0; i < cols; i++) {
        eigenvectors[i] = (double*)malloc(cols * sizeof(double));
    }

    // Data processing steps
    center_data(data, rows, cols, mean);
    covariance_matrix(data, rows, cols, cov);
    jacobi_eigen(cov, eigenvalues, eigenvectors, cols);

    // Calculate total variance
    double total_variance = 0.0;
    for (int i = 0; i < cols; i++) {
        total_variance += eigenvalues[i];
    }

    // Sort eigenvalues and eigenvectors in descending order
    for (int i = 0; i < cols - 1; i++) {
        for (int j = i + 1; j < cols; j++) {
            if (eigenvalues[i] < eigenvalues[j]) {
                // Swap eigenvalues
                double tmp = eigenvalues[i];
                eigenvalues[i] = eigenvalues[j];
                eigenvalues[j] = tmp;

                // Swap eigenvectors
                for (int k = 0; k < cols; k++) {
                    double t = eigenvectors[i][k];
                    eigenvectors[i][k] = eigenvectors[j][k];
                    eigenvectors[j][k] = t;
                }
            }
        }
    }

    // Calculate explained variance ratio
    double explained_ratio1 = eigenvalues[0] / total_variance * 100;
    double explained_ratio2 = (eigenvalues[0] + eigenvalues[1]) / total_variance * 100;

    // Output results
    printf("\n=== Principal Component Eigenvalues ===\n");
    for (int i = 0; i < cols; i++) {
        printf("PC%d: %.5f (%.2f%% of variance)\n",
               i + 1,
               eigenvalues[i],
               eigenvalues[i] / total_variance * 100);
    }

    printf("\n=== Explained Variance ===\n");
    printf("First principal component: %.2f%%\n", explained_ratio1);
    printf("First two principal components: %.2f%%\n", explained_ratio2);

    printf("\n=== PCA Results (First Two Principal Components) ===\n");
    printf("%-8s %-12s %-12s\n", "Sample", "PC1", "PC2");
    printf("----------------------------------------\n");
    for (int i = 0; i < rows; i++) {
        double pc1 = 0, pc2 = 0;
        for (int j = 0; j < cols; j++) {
            pc1 += data[i][j] * eigenvectors[0][j];
            pc2 += data[i][j] * eigenvectors[1][j];
        }
        printf("%-8d %-12.4f %-12.4f\n", i + 1, pc1, pc2);
    }

    // Free allocated memory
    free(mean);
    for (int i = 0; i < cols; i++) free(cov[i]);
    free(cov);
    free(eigenvalues);
    for (int i = 0; i < cols; i++) free(eigenvectors[i]);
    free(eigenvectors);
}

int main() {
    int rows, cols;
    double **data;

    // Please modify the path to your CSV file
    const char* filename = "yield_df.csv";
    printf("Attempting to open file: %s\n", filename);

    data = read_csv_numeric(filename, &rows, &cols);
    if (data != NULL) {
        printf("Successfully read %d rows, %d numeric columns.\n", rows, cols);
        pca(data, rows, cols);

        // Free allocated data memory
        for (int i = 0; i < rows; i++) free(data[i]);
        free(data);
    } else {
        printf("Failed to read CSV file.\n");
    }

    return 0;
}
